function f = covar(theta2,resid,sj,v,x2,p,QI,ind_market,A1,W,d_d,ztilde,xtilde,deltaold,sales_tax)

% this function computes the var-covar matrix; 

N=length(resid)/2;
k2=length(theta2); % number of nonlinear parameters
z_D=ztilde(1:N,1:d_d); % demand instruments 
z_S=ztilde(N+1:end,d_d+1:end); % supply instruments 
d_S=size(z_S,2);
xi=resid(1:N,1);
omega=resid(N+1:end,1);


h=1e-8;
H=eye(k2,k2)*h;
ai=computey(theta2+H(1,:),sj,v,x2,p,QI,ind_market,A1,deltaold,sales_tax);
af=computey(theta2-H(1,:),sj,v,x2,p,QI,ind_market,A1,deltaold,sales_tax); 
K=(ai-af)/(2*h); 
if k2>1
    for i=2:k2
        Ai=computey(theta2+H(i,:),sj,v,x2,p,QI,ind_market,A1,deltaold,sales_tax);
        Af=computey(theta2-H(i,:),sj,v,x2,p,QI,ind_market,A1,deltaold,sales_tax);
        K=[K (Ai-Af)/(2*h)];
    end
end

Gamma=(ztilde'*[-xtilde K]);

f=inv(Gamma'*W*Gamma); 